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ABSTRACT 

We use numerical hydro dynamic simulations to investigate prestellar core for- 
mation in the dynamic environment of giant molecular clouds (GMCs), focusing 
on planar post-shock layers produced by colliding turbulent flows. A key goal is 
to test how core evolution and properties depend on the velocity dispersion in 
the parent cloud; our simulation suite consists of 180 models with inflow Mach 
numbers Ai = v/cs = 1.1 — 9. At all Mach numbers, our models show that 
turbulence and self-gravity collect gas within post-shock regions into filaments 
at the same time as overdense areas within these filaments condense into cores. 
This morphology, together with the subsonic velocities we find inside cores, is 
similar to observations. We extend previous results showing that core collapse 
develops in an "outside-in" manner, with density and velocity approaching the 
Larson-Penston asymptotic solution. The time for the first core to collapse de- 
pends on Mach number as tcou oc A^^^^^pg for po the mean pre-shock density, 
consistent with analytic estimates. Core building takes 10 times as long as core 
collapse, which lasts a few x 10^ yrs, consistent with observed prestellar core life- 
times. Core shapes change from oblate to prolate as they evolve. To define cores, 
we use isosurfaces of the gravitational potential. We compare to cores defined 
using the potential computed from projected surface density, finding good agree- 
ment for core masses and sizes; this offers a new way to identify cores in observed 
maps. Cores with masses varying by three orders of magnitude (~ 0.05 — 50Mq) 
are identified in our high-Al simulations, with a much smaller mass range for 
models having low Al. We halt each simulation when the first core collapses; at 
that point, only the more massive cores in each model are gravitationally bound, 
with Eth + Eg < 0. Stability analysis of post-shock layers predicts that the first 
core to collapse will have mass M oc w ~^/^Pq ^^^T^/^, and that the minimum mass 
for cores formed at late times will have M oc p^^^'^T'^ , for T the temperature. 
From our simulations, the median mass lies between these two relations. At the 
time we halt the simulations, the M vs. v relation is shallower for bound cores 
than unbound cores; with further evolution the small cores may evolve to become 
bound, steeping the M vs. v relation. 
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Subject headings: ISM: clouds — ISM: globules — stars: formation 



Introduction 



Star formation begins with the creation of dense molecular cores, and und erstanding 



how cores grow and evolve is es s ential to identifyin g the origin of stellar properties (IShu et al. 



19871 : iMcKee fc Ostrikerl 120071 : lAndre et al.ll2008f ). Through the 1990s, the prevailing the- 
oretical picture was of slow core formation and evolution mediated by ambipolar diffusion, 
followed by core collapse initiated from a qua sistatic, centrally-concentrated state (e.g., 
Mouschoviad Il987l : iMouschovias fc Ciolekl Il999l ). Current observations, however, indicate 
that r nagnetic field strengths are insufficient to provide the dominant support of molecular 
cores (ITroland fc Crutcheii 120081 ). In addition, over the past decade, a conception of star 
formation has emerged in which supersonic turbulence drives struc ture and evolution within 



giant molecular clouds (GMC s) on a wide range of scales (e.g., iBallesteros-Paredes et al. 



20071 : iMcKee fc Ostrikerll2007| ). Because supersonic turbulence can compress gas to densities 
at which gravitational collapse can rapidly occur, it is likely to be important in the initiation 
of prestellar cores. Ultimately, models of core formation and evolution mu st take into ac- 
count both moderate magne tic fields (with diffusion) and strong turbulence (IKudoh &: Basu 
20081 : iNakamura &: Lill2008l ). In order to gain insight into the physics involved, however, it 
is informative t o focus on individual lim iting cases and explore dependence on parameters. 
Here, following iGong fc Ostriken ( l2009l ) but generalizing to three dimensions, we consider 
core building and evolution in the turbulence-dominated, unmagnetized limit. 

Observations of dense cores in GMCs have provided detailed information on indi- 
vidual core properties as well as statistics of core populations (see e . g., the reviews of 



di Francesco et al.l 120071 : IWard- Thompson et al.l 120071 : iBergin fc Tafallal 120071 : I Andre et al. 



20081 ). These properties, including internal structure and kinematics, durations of differ- 
ent evolutionary stages, and distribution of core masses, constrain core formation theories. 
In terms of structure, cores are observed to be centrally concentrated at all stages, with 
the specific profile fits differing depending on the stage of evolution. Cores can gener- 
ally be fit w'ith a uniform- density inner region surrounded by a power law o c (e.g.. 



Shirlev et al.ll2000l : iBacmann et al.ll2000l : lAlves et al.ll2001l : iKandori et al.l 120051 : iKirk et al. 
2OO5I ): this shape is consistent wit h expectations for both sta tic Bonnor-Ebert (BE) pressure- 



supporte d isothermal equilibria (Bonnor 19561 EbertI 19551), and for collapsing isothermal 
spheres ( Bodenheimer fc Sweigartl 19681 : Larso 2 Il969[ iPenstonI Il969l ) . The center-to-edge 
density contrast is frequently larger than the maximum possible for a stable BE sphere, 
however, and the inferred temperatures based on static BE fits are also often larger than 
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observed te mperatures. Although in prin ciple some support could be provided by magnetic 



fields (e.g., ICiolek fc MouschoviasI Il994l ). ano ther possibility is that these "supercritica l" 



cores are in fact collapsing rather than static ( iDapp &: Basull2009l : iGong fc Ostrikerll2009[ ) 



In terms of kinematics, dense, low- mass cores generally have subsonic in ternal veloc- 



Goodman et al. 


1998; 


Caselli et al. 


2002 


2007; 


Lada et al 


. boos 


). Some prestellar 



Mverd Il983l: 



Andre et al. 



tions throughout their interiors based on asymmetry of rn olecular lines that trace dense gas 



[e.g., 



Lee fc MyersI 119991 : iLee et al.ll200ll ; ISohn et al.l 120071 ). For cores containing protostars, 
signa tures of supersonic inward motions on small scal es (~ 0.01 — 0.1 pc) have been observed 



Gregersen et al. 1997; Di Francesco et al. 200l[): these are believe d to be indicative of 



(e.g., 

gravitationally-induced infall. In very recent work, iPineda et al.l (|2010( ) have used NH3 ob- 
servations to identify a sharp transition from supersonic to subsonic velocity dispersion from 
outer to inner regions in the core B5 in Perseus. 

Several recent statistical studies have reach ed similar conclusions regard i ng the dura 



tions of successive stage s of core evolution (e.g., IWard-Thompson et al.l 120071 : Enoch et al. 



2008 



Evans et al.ll2009l ). with prestellar and protostellar (class 0) stages having compara- 
ble lifetimes. The typical duration for each of these stages is a few times the gravitational 
free-fall time 

Q \ 1/2 



4.3 X 10° yr 



10^ cm- 



-1/2 



at the mean core density p = lAmfjnH, amounting to ~ 1 - 5 x 10^ yr for typical conditions. 
With prestellar lifetime s considerably below the amb ipolar diffusion time for strong magnetic 



field tAD ~ 10^// (s-g- Mouschovias fc Ciolek 19991), t his suggests that observed cores are 



trans-critical or supercritical (see Ciolek fc Basu 2001) with respect to the magnetic fieldjll 
This conclusion is also supported by magnetic field Zeeman observations ( jTroland fc Crutcher 
20081 ). indicating that cores have mean mass-to-magnetic- fiux ratios twice the critical value. 
Thus, magnetic field effects appear to be sub-dominant in terms of supporting cores against 
collapse, and ambipolar diffusion does not appear to control the dynamics of core forma- 
tion and evolution. As magnetic fields are non- negligible, however, magnetohydrodynamic 
(MHD) stresses may still affect GMC and core dynamics. 



Empirical measurements of core mass functions (CMFs) (e.g.. 




otte et al.lll998l: 


Testi & Sargent 


1998; 


Johnstone et al. 


2000, 


2001; 


Motte et al. 


2001; 


Onishi et al. 


2OO2I: iBeuther & Schilkel 



1 The critical mass-to-magnetic-fiux defines the minimum that perm i ts Kravitational collapse in t he field- 
freezing limit (e.g. Mestel fc Spitzer 19561 : Mouschovias fc Spitzei 19761 : Nakano fc Nakamura 19781 ). 
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2004; 


Reid & W] 


son 


2005 




2006; 


Stanke et al. 


2006; 


Enoch et al. 


2006 




4.1ves et al. 


2007 




Ikeda et al.lboO?. 


2009 


: Ikeda & Kitamurall2009: 


Nutter & Ward-Thompson 


2007; 


Simpson et al. 


2008|; 


Konvves et al. 2 


010) show that CMFs have a remarkable similarity in shape to stellar 



initial mass functions (IMFs, see e.g. iKroupa 



lower mass by a factor of 3 - 4 (see e.g., lAlves et al. 



2001 



2007 



Chabrier 20051) . with a shift toward 



Rathborne et all 120091 ). The 



characteristic/turnover mass of observed CMFs ranges from 0.1-3 M©, although there are 
uncertainties in this associated with lack of spatial resolution at the low mass end. 

Many theoretical eff orts have cont ribut ed to i r iterpr eting the observed properties of 
cores. The classic work of iBonnorl (jl956[ ) and lEbertl ( 119551 ) provided the foundation of later 
studies, by determining the maximum mass of a static isothermal sphere that is dynamically 
stable. In terms of the boundary pressure Pedge = Pedgecf or mean internal density p = 
2.5pedge) this maximum stable mass is 

-1/2 / \ 3/2 

lOKJ ■ 



M, 



BE 



1.2- 



1.9- 



2.3 Mcr 



(G^Pedgc)^/' ^■^(G3p)l/2 ^'^® Vl04cm-3 

{kT/ fiY^'^ is the internal sound speed in the core. 



Here, c. 



(2) 



1969; 



2003 



of individua 


, pre-existing cores 


Bodenheimer & Sweieart 1968; . 


jarson 1969; 


^enston 


Hunter 1977 


Foster & Chevalier 


1993 


Odno et al. 


1999 


Hem 


aebe 


le et al. 


2003 


Motovama & Yoshid 


Vorobvov & Basu 2005: Gomez et al. 


2007; 


Burkert & AlvesI 


200( 


}). These simulations 



include initiation from static configurations that are unstable, and initiation from static, 
stable configurations that are subjected to imposed compression, either from enhanced ex- 
ternal pressure or a converging velocity field, or a core-core collision. A common feature of 
the results is that the collapse generally starts from outside and propagates in as the cen- 
tral density increases. At the time of singularity formation, the density profile approaches 
the "Larson-Penston" asymptotic solution p = S.SQc^/iAir Gr'^) and the c entra l velocity is 
comparable to the value — 3.28cs derived by iLarsoru (119691 ) and iPenstonI (119691 ). However, 
these previous studies have not considered core evolution within the larger context, in par- 
ticular including the process of core formation. Since the formation process may affect later 
evolution, it is important to develop unified models. 

At GMC scales, a number of groups have investigated the CMFs that result from numer 



ical simulations of 


turbulent, self-gravitating systems (see e.g., 


Klessen 


2001 


Gammie et al. 




2003 


Bonnell et a 


2003; 


Li et al. 


2004 


Tillev & Pudrit2 


! 2004: 


Heitsch et al.l 


2008; 


Clark et al. 


2008 


loffner et al. 


20081; 


Basu et al. 


2009: Smith et al. 


2009). 


These models have shown - 



for certain parts of parameter space - features that are in accord with observed CMFs: mass 
functions dominated by the low end with a peak and turnover near 1 Mq, and a high-mass 
power-law slope (at least marginally) consistent with the Salpeter value. These simulations 
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have not, however, had sufficient resolution to investigate the internal properties of indi- 
vidual cores that form. In addition, these studies have not quantified how the core masses 
depend on the large-scale properties of the turbulent medium (see below). 



Taking the previous numerical simulations of individual cores one step further. iGong fc Ostriker 



( I2OO9I ) initiated a study of dynamically induced core formation and evolution in supersonic 
converging flows, focusing on the spherical case. In these simulations, the density is initially 
uniform everywhere: no initial core structure is assumed. Instead, dense cores form inside 
a spherical shock that propagates outward within the converging flow. Over time, cores 
become increasingly stratified as their masses grow. Eventually, the core collapses to create 
a protostar following the same "outside-in" pattern as in models initiated from static con- 
dition s. Subsequently, the den se envelope falls into the center via an inside-out rarefaction 
wave ( IShij|l977l : iHunterl 119771 ): this is followed by a stage of late accretion if the converging 
flo w on large scales c ontinu es to be maintained. The unified formation and evolution model 
of iGong fc Ostrikerl (120091 ) explains many observed core properties, including BE-sphere- 
like density proffies, subsonic internal velocities w ithin cores, and short co re lifetimes with 
comparable prestellar and protostellar durations. iGong fc Ostrikerl (120091 ) also found that 
the inflow velocity of the converging flows affects core lifetimes, masses, sizes and accretion 
histories. Realistic supersonic inflows in clouds are not spherical, however, while mass inflow 
rates are affected by geometry. Thus, the quantitative results for masses, lifetimes, etc., as 
a function of Mach number and ambient density may differ for more realistic geometry. 

Numerical results on core formation have not reached consensus on how the character- 
istic mass in the CMF, Mc, depends on the bulk properties of the cloud - its mean density 
Po = (p), sound speed c^, and turbulent velocity dispersion fturb- Some have suggested that 

in the CMF (e.g,. ' 



the Jeans mass of the cloud at its mean density (M j = c^7r^/^(G^po) ^^^) determines 



Klessen 



2001 



Bonnell et al.ll2006f). whil e others have found values of M, 



well b elow Mj (see e.g.. lGammie et al.ll2003l : iLi et al.ll2004j ). As noted by lMcKee fc Ostrikei 



( 120071 ). the difference between these conclusions is likely related to the Mach number of 
turbulence: the value found for M^/Mj is lower in simu lations wher e the Mach number 
M. = Vtnrh/ Cs IS higher. Indeed, more recent simulations by I Clark et al.l (|2008[ ) provide some 
indication that increasing Ai lowers the value of in the CMF; they did not, however, 
conduct a full parameter study. 

Supersonic turbulence makes the density in a GMC highly non-uniform, creating a 
log-normal probability distribution function (PDF) in which most of the volume is at den- 
sities below pn and most of the mass is at densities above po (e.g.. IVazquez-Semadenilll994j : 
Padoan et al.lll997l : lOstriker et al.lll999[ ). Given that the l og-normal PDF allows for a range 
of Jeans masses (or Bonnor-Ebert masses; Mbe oc Mj), iPadoan fc Nordlundl (120021 12004[ ) 
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proposed that t he CMF is set by div iding: the total available gas mass at each density into 



unstable cores. iPadoan et al.l ( 120071 ) propose that the peak mass in the CMF is given by 
Mc = 3Mbe,o/M^^ for Ma = Vmrh/vA the Alfven Mach number in a cloud, and Mbe,o 
the Bonnor-Ebert mass evaluated at the mean cloud density uq. Here, va = B / [AnpY/"^ 
is the Alfven speed. For realistic mean GMC density no ~ 100 cm~^ and M.a ~ 1 — 4, 
from Equation ([2]) the Padoan et al formula in fact yields > 15 M©; only if one chooses 



case, 



Padoan et all (120071 ) propose that = AMbe,o/M 



1.7 



Hennebelle & Chabriei 


( 


2008[) 


Lgnetized case bv Padoan et al. 



( 20071 ). and advocate a formula similar to their unmagnetized one: Mc ~ Mbe,o/-M^^'^- Since 
Ai > 10 in massive GM Cs, the se formulae yield more r ealist ic values Mc ~ Mq. Neither 



the iPadoan et al.l (120071 ) or the iHennebelle &: Chabrierl (120081 ) proposal has, however, been 
tested directly using self-gravitating numerical simulations. 

In this contribution, we present results on core formation and evolution based on a large 
suite of 3-dimensional numerical simulations. Each simulation models a localized region of 
a turbulent cloud in which there is an overall convergence in the velocity field. Under the 
assumption that there is a dominant convergence direction locally, we choose inflow along a 
single axis, so that convergence is planar. With the more realistic geometry afforded by the 
current simulations, we are able to check the results obtained by iGong &: Ostrikerl (120091 ) 
for core building and collapse in supersonic flows. We are also able to explore how the 
characteristic core mass is related to the velocity of the converging flows. Since the speed of 
converging flow is assumed to reflect the amplitude of the largest-scale (dominant) motions in 
a GMC, this relates the characteristic core mass to the turbulent Mach number in its parent 
GMC. Although a number of previous studies of core formation have been conducted, the 
present investigation is distinguished by our systematic study of Mach number dependence, 
together with our focus on internal structure and kinematics of the cores that form. 

The plan of this paper is as follows: In Section 2 we provide a physical discussion of 
self-gravitating core formation in the post-shock dense layers, identifying the mass, size, 
and time scales expected to be important. In Section 3, we summarize the governing 
equations and methods used in our numerical simulations. Section 4 describes the devel- 
opment of core structure and evolution in our models, paying particular attention to the 
influ ence of Mach number jy i on the evolution, and comparing collapse of individual cores 
with iGong fc Ostrikerl (120091 ) . Section 5 describes our method of core- finding, in which the 
largest closed contour of the gravitational potential determines the core size. We demon- 
strate that this method can be used for both three dimensional and two dimensional data 
with similar results, and can thus be applied to find cores in observed clouds. Section 6 
describes the relations between core properties (core mass, core radius and core collapse 
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time) and the large-scale Mach number of the converging flow, relating to the expectations 
from gravitational instability discussed in Section 2. In Section 6, we also quantify core 
shapes, and explore the relationship between core structure and kinematics. Section 7 sum- 
marizes our new results and discusses our findings in the context of previous theories and 
observations. 



2. The characteristic core mass and size 

Prior to describing our numerical model prescription and results, it is useful to sum- 
marize the scales that are likely to be relevant for formation of self-gravitating cores in 
GM Cs. We shall assu me approximately isothermal conditions, consistent with observations 



[e.g. 



Blitz et al.l 120071 ). The isothermal sound speed at a temperature T is 



y \ 1/2 



c, = 0.20 km s-' l^-^j . (3) 

If the density within clouds were uniform, the spatial scale relevant for gravitational insta- 
bility would be the Jeans length 

^Gio) ='-''pnTow^J ItokJ ' 

evaluated at the mean density po- The corresponding Jeans mass is 

Note that Pq{Lj/2)^ or po4vr(Lj/2)^/3 is sometimes used for the Jeans mass. The Bonnor- 
Ebert mass (eq. [2]) for Pedge = -Po = Poc^ is Mbe = 0.22Mj(po)- The Jeans time at the mean 
cloud density is 

We shall use the Jeans length, mass, and time at the unperturbed density as our code units 
of length, mass, and time: Lo = Lj, Mq = Mj, and to = tj. 

Of course, GMCs are highly inhomogeneous, with core formation taking place in the 
overdense regions that have the shortest gravitational times. If the overdense regions within 
GMCs are produced by shocks in the turbulent, supersonic flow, their density, and therefore 
the mass scale and length scale for growth of self-gravitating structures, will be related to the 
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shock strength. Strongly magnetized shocks have less compression than weakly magnetized 
shocks (while both will be present in a turbulent flow), so we concentrate on the latter case. 

If gravitationally unstable cores develop only in gas that has been strongly compressed 
by shocks, the actual bounding pressure will be much larger than Pq = pocl- In particular, 
an isothermal shock with Mach number Ai will produce a post-shock region with pressure 
-Ppost-shock = Po^^ = -^^PoCs ^ Pq. Thus, if cores preferentially form in stagnation regions 
between shocks of Mach number Ai, then one can define an effective Bonnor-Ebert mass for 
these core-forming regions within the turbulent flow by setting Pedgc = -Ppost-shock in equation 

MB^,po.,_3hock = ^-^jGW^Ji = ^® (lk£^) (lOW^) ^ (fk) • 

The above simple argument suggests M oc v^^p^^^'^T'^ for the minimum mass of a star that 
forms via collapse of a core in a turbulent cloud with velocity dispersion v, mean density po, 
and temperature T. 

Equation ([7]) provides a mass scale for fragmentation within post-shock regions, but 
in fact instabilities take some time to develop. Thus, it is useful to consider the evolution 





31meereen & Elmeereen 


1978 


1994 


. Iwasaki & Tsuribe 


200^ 





9781 : iLubow fc Pringlelll993l : IVishniadll994l : IWhitworth et al. 



For inflow Mach number Ai, the surface density of the post-shock layer at time t is 

^it)= poiv,,+ -v,,.)t = 2poMcst, (8) 

where Vz,+ and Vz~ are the upward and downward converging velocities. If the sheet is 
not vertically self-gravitating, its half-thickness is if = S(t)/2pp where pp ~ PoA^^ is the 
post-shock density. The non-self-gravitating half-thickness is thus 

2po M Cs t _ cj^ 
"^^"^ 2po-M2 ~ M- ^ ' 

As the surface density of the sheet increases, self-gravity will become increasingly important 
in confining the gas. In the limit of hydrostatic equilibrium, the height approaches 

c 

^ ^ ^ 2710 p^Mf 

Note that the transition from non-self-gravitating (if nsg oc t) to self-gravitating (iJgg oc t^^) 
occurs at a time near 
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defined by the condition H^g = Hj 



nsg- 



The disper sion relation for in-plane modes in a slab, allowing for non-zero H (e.g. 



Kim et al.ll2002[ ). is 



2 2 7 2 lirGYjk /-, 



For the critical mode = 0, so that 



where 



k,,,,H (1 + fcerit^) = 27rif ^ = 27r-^, (13) 



^J,2D = ^ (14) 



is the Jeans length for an infinitesimally-thin layer. The solution to equation (fT3|) is 

_ 27r 2 _ A-KGpQtM 2 

Kcrit - T - -T7^ - - - -T7^, 



(15) 



'J, 2D J \ ij,2D 

so that 

/ X 1/2 / X 1/2 

1+ l + 87ryi^ ^ 1+ I + Stt-^ 



-^"'^ - 2 2GMM ^ ■ ^^^^ 



The corresponding critical mass (Aci.it/2) S is 



M - 



l + (l + 87rT^)i/2l' 



^ crit 



(17) 

32G2poA^ t 

Note that H/ Lj_2D initially increases in time, during the non-self-gravitating stage {H^sg/ -^j,2D 
2Gpot'^), and then approaches a constant (ifsg/-Z^j,2D = 1/^)- any time, all wavelengths 
A > Acrit have < 0, so that overdense regions of the corresponding sizes and masses 
M > Merit grow relative to their surroundings. 

During the non-self-gravitating stage, the critical mass has a minimum value at time 

3 N^/^ 



given by 



'^crit,nsK,min 1 i^ ^ I 0.14tj O.Gltgcr 

\l67rGpo ' 



_ 3V3^ 1 

Mcrit,nsg,min - g (G^p^y/^M ^^^^ 
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The numerical coefficient in equation ( !T9|) is 1.15; note that this is almost the same as in 
equation ([7]). 

At late time, the critical mass from Equation f[T7|) with -ffsg/-^j,2D = 1/vr becomes 

^"^''^^ " 2G^PoMt ^ ^^^^ 

Expressing Mcrit,sg in terms of the virial parameter a^n = 5(t^-R/GMgmc of the GMC, and 
using (Jt, = Aics and Mgmc = tt^R^Sgmc = 47r_R^po/3, we have 



Mcnt,sg- 1^ 20 j (G^poY/'M E ■ ^ ^ 

Here ay is the large-scale one- dimensional velocity dispersion in GMCs, which will be re- 
sponsible for the largest scale, strongest shocks. Taking Ovir = 2, the coefficient in Equation 
fl22|) is 0.97, so this is very similar to equations ([7]) and f lT9|l if S ~ Sgmc- In dimensional 
units, the critical mass (for S = Sgmc) is 



M„itsg = 2.5 Mq "' J • (23) 

cnt,sg oVlkms-i/ V102cm-3y VlOK/ ^ ^ 

As noted above, equations ([7]), f fT9|) - (20) and fl22|) - fl23|) all have a similar form. An 
important task for numerical simulations is therefore to test the hypothesis that the char- 
acteristic mass scale of collapsing cores formed in turbulent, self-gravitating GMCs follows 
this scaling, i.e. 

where is a dimensionless coefficient. 

The critical mass given above is the smallest mass that can collapse, given infinite time. 
Since the growth rate depends on scale (and is formally zero for critical perturbations), at 
any finite time only cores that have grown sufficiently rapidly will be nonlinear enough to 
collapse. It is therefore useful to consider how much growth has occurred at a given time. 
Consider a perturbation of wavenumber k that instantaneously has d^5E/dt^ = — t<j^5S so 
that 5S = 5Einite^ where T = In (5S/5Sinit) = Ji-u^Y^'^ dt. Using equation 
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where tmin is the instant when S is large enough that perturbations of wavenumber k start 
to grow (— ^ 0). With S = 2pQCsMt, tmin = Csk{l + kH) / {AnG PqM) . If we assume 
kH ^ 1 (see below), then 

V='^-^n''\r-n/2f'\ (26) 
where k = kcg/ \/2tiGpq7A and r = t^/^jxGpoM.. 

At a given time t (or r) during the evolution, the mode k^ (or /«„) that has grown the 
most has d\nT/dk = 0, which gives 

f^m = (27) 

and Fjnax = r(^m) = y/^i^m = V^ttG poAdt'^ /2. The mass of this most- amplified mode is: 



where the time is 



2r , \ / T X 1/2 



and fc„ = {r^.^J^/3)^/'^{27rGpoMY/'^ /cs, so that 

1/2 



With Fmax = 1, the numerical coefficient for in Equation ( l28|) is 3.30, and Equation 
f l2I?]) gives t = 0.34tjA^~i/2, corresponding to r = 1.5. Note that for low Mach number, this 
time exceeds tgg (see eq. [TT]) . whereas for high Mach number it does not. Also, note that 
with H < CstJM = H^^^ (see eqs. [9] -[11]), k„,H < k^H^^^ = tUL{VSM)-^/^ . Taking 
Tmax = 1, kmH < 0.8 for > 1, with kmH <^ 1 for ^ 1. This verifies self-consistency 
of the assumption made in obtaining equation (1261) . 



Written in terms of v,po, and T, the most-amplified mass is 

M. = 19.1 Me f^^)"'' ( .V ( -^Y' (r^nax)-^/^ ■ (31) 

°Vlkms-i/ Vl02cm-3y \10KJ ^ ^ ^ ' 

Comparing equation ( 1ST]) with equation ( 12^ . we see that a different dependence on velocity 
(or Mach number) is expected for the first core to collapse (equation EI]), compared to 
the typical core to form eventually feauation [231). Similar res u lts to equation f[28]) have 
previously been discussed by other authors. [Whitworth et al.[ (Il994j ) point out that the 
fastest-growing scale ~ -Z^j,2D ~ c^/ {Gp^M-t) will become nonlinear if the time exceeds the 
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growth time ~ Lj^2d/cs ~ {GpoJ^t)~^, which occurs for t ~ (GpoAi)^^^'^ (cf. our eq. 
[29]) . This corresponds to a length scale ^fragment ~ Cs{GpoAi)~^^'^ (cf. our eq. [30]), and 
a mass scale Mfragment ~ c^{G^poAi)~^^'^ (cf. our eg. 1251). By direct i n tegrat ion of the 
perturbation equation of the converging-flow system, llwasaki fc Tsuribd ( 120081 ) find that 
the fastest-growing mode becomes nonlinear at time .965^^'^ {G poAi)~^^^ , for the initial 
amplitude (cf. our eq. [291 which has a coefficient 0.6 if Fmax = !)• 



Finally, we note that the characteristic mass scale at late times given in equation ( I24p can 
be connected to observed core mass scales using the empirical relationships among turbulence 
level, size, and mass for GMCs. In terms of the viral parameter a^v = 5a^R/ (GMgmc) and 
the GMC surface density Sgmc = 4po-R/3, equation f[M|) can be re-expressed as 

Sgmc 



Mc = 1.5^ 



X 1 M0 



10^ j 



100 Mq pc- 



1/2 



(32) 



With ftvir ^ 



1 - 2 and EpMc ~ 100 Mfn pc for observed clouds ([Solomon et al.l [1987 



McKee fc Ostriker[ [20071 : [Heyer et al.l [2009[ ). the mass scale is intriguingly similar to the 
characteristic (peak) mass of CMFs within nearby molecular clouds. This relation potentially 
also offers a prediction for the peak of the CMF (and ultimately the IMF) when stars 
form under conditions different from those in most Milky Way GMCs. In particular, high 
temperature (up to ~ 70 K) may hold in starburst regions where the radiation field is strong 
and turbulent dissipation rates are high; since the temperature dependence of equation ([32]) 
is steeper than the dependence on surface density, this could imply higher masses under 
those conditions. 



3. Methods for numerical simulations 



The numerical simulations we present here are cond ucted with the Athena co de ([Gardiner &: Stone 
20051 . 12OO8I : Istone et allbood : Istone fc Gardinei]l2009h . using the HLLC solver jTorolll999h 
and second order reconstruction (IStone et al.|[2008l ). To calculate the self-gravity of our slab 
domains, which are periodic in-plane an d open in the z dire c tion, the Fast Fourier Trans- 
formation (FFT) method developed by [Koyama fc Ostriker[ (|2009[ ) is used. We solve the 
three-dimensional equations of hydrodynamics. 



dp 
dt 



(33) 



(34) 
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and the Poisson equation, 

V2<|> = AttGp, (35) 
where $ is the gravit ational potential The isothermal assumption P = cj.p is adopted. 



Pavlovski et al.l ( 120061 ) found the isothermal approximation is adequate for simulations of 



the interstellar medium even with strong turbulence, which implies strong shocks in GMCs. 

The code unit of density po is a fiducial density representing the volume-averaged am- 
bient density in a cloud on large scales; this characterizes the mean density of converging 
flows. For the code unit of velocity, we adopt the isothermal sound speed Cg (see eq. [3]). For 
the unit of length, we adopt Lq = Lj, the Jeans length at the fiducial density (see eq. H]). 
The mass and time units for the simulation are then Mq = Mj (see eq. [5]) and to = tj (see 
eq. E]). 

In making comparison to observations, the total surface density integrated through the 
domain 

p dz 



p{x,y,z)dz = T.o —— (36) 
J Po 

is useful, for Sq = poLj = 9.49 M© pc-^{T/10KY/^{nH,o/10^ cm-^y/'^. In terms of the 
column density of hydrogen, 

iV„ = ^ = A^„/^^ (37) 

for A^o = uqLj = 8.51 x 10^° cm-2(r/10ir)i/2(n//,o/10^ cm-3)V2. The mean line-of-sight 
velocity is calculated by 

, , / pviosds 

and the corresponding dispersion of (fios) is defined as 

2 J p{vios - {vios)yds 
j pas 

where ds = secO dz and 9 is the tilt angle of the observer with respect to the z axis. 

Our model prescription consists of a converging flow augmented with turbulent velocity 
perturbations. In our parameter survey, the Mach number M. of the inflow velocity ranges 
from 1.1 to 9. Thus, two flows converge toward the central plane 2; = from the upper z- 
boundary (with mean velocity —Aicg) and the lower 2;-boundary (with mean velocity M.Cs)- 
The initial density is uniform and set to po; and the density at the inflowing ^-boundaries 
is also set to po throughout the simulation. The boundaries in the x and y directions are 
periodic. 
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For both the whole domain initially and the inflowing gas subsequently, we apply per- 
turbations following a Gaussian random distribution, with a Fourier power spectrum of the 
form 

{\6v,\') oc k-\ (40) 

for \kL/27r\ < N/2, where is the resolution and L is the size of the simulation box 
in X an d y. The power spectru m is appropriate for supersonic turbulence as observed in 



GMCs (IMcKee fc Ostrikerl 120071 ) . The perturbation velocity fields are pre-generated with 
resolution 256^ in a box of size L^. The perturbation fields are advected inward from the 
2;-boundaries at inflow speed M-Cs- at time intervals At = /S.z / {M.Cs), slices of the pre- 
generated perturbation fields for Vx, Vy and Vz are read in to update values in the ghost zones 
at the z-boundaries. 

In addition to exploring dependence on the mean inflow Mach number M., we also 
test dependence on the ampli tude of turbulent perturbations on top of this converging flow. 
From the scaling law (see e.g.,|Larso nl ll98l[lHever fc Bruntl[20oi ) of self-gravitating molecular 



clouds, 5v{l) oc we can write the velocity dispersion at scale / in terms of cloud-scale 
one- dimensional velocity dispersion and cloud radius R as 5vi£,{l) = c^^,(//2i?)^/^. The 
velocity dispersion at the scale of the simulation box L is 

Cs c, \2Rj c, [LjJ [LjJ ■ ^ ^ 

In terms of the viral parameter a^ir = 5 cr^R/ (GM), where M = AttR'^Pq/S is the cloud mass, 
the ratio between ay and Cs is 



Cs V 15 / Lj 

Solving equation P2|) for 2R/Lj and substituting into equation (HTj) . we have the amplitude 
of perturbation for the simulation box: 



15 J UJ \Lj 



(43) 



Thus, if the size of the simulation box is L = Lj and a^r = 1^2, the perturbation amplitude 
would be 

Cs \CsJ 

If we take the Mach number of the inflow, A4, as comparable to the value (Jy/cs of the whole 
cloud, then equation (jUj) implies that higher converging velocities would be associated with 
higher amplitudes for the perturbation fields, for a given simulation box size Lj. To test 
the influence of the perturbation amplitude, we conduct two sets of simulations with 10% 
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and 100% of the value 5vid{Lj) = (Ai/SY^'^Cs. Hereafter, we denote these cases as low 
amplitude and high amplitude initial perturbations, respectively. 

For each Mach number M at each amplitude, we run 20 simulations with different 
random realizations of the same perturbation power spectrum, in order to collect sufficient 
statistical information on the core properties that result. The whole set of simulations there- 
fore consists of 180 separate runs. The resolution for low amplitude perturbation simulations 
is Nj. X Ny X Nz = 256 x 256 x 96, with domain size x Ly x Lz/Li^j = 1 x 1 x 0.375; 
for high amplitude the resolution is x Ny x = 256 x 256 x 160, with domain size 
LxX LyX Lz/Li^j = 1 X 1 X 0.625. The domain in the z direction is smaller than in the x and 
y directions since the reversed shock generated by the inflow only propagates a relatively 
short distance and the post-shock dense layer is thin, i.e., the basic geometry remains planar. 
The domain in the z direction is large enough so that the post-shock layer does not evolve 
to reach the z boundaries. 

We note that our assumption of perturbed velocities but uniform densities in the in- 
flowing gas is not fully realistic, since the flow entering a strong shock within a GMC will 
in general have internal density structure@ In fact, the velocity perturbations we introduce 
do lead to moderate (order-unity) density fluctuations, as we have found by conducting 
comparison simulations with self-gravity turned off. These density fluctations are what seed 
the growth of self-gravitating structures. The main emphasis of the current work is to in- 
vestigate how the development of self-gravitating structures depends on the inflow Mach 
number, which sets the mean density (and hence the gravitational timescale) in the post- 
shock layer; previous studies have not tested the Mach number dependence of gravitational 
fragmentation. By varying the velocity perturbation amplitudes of the inflow, we have begun 
to explore the effect of pre-existing density structure on self-gravitating core development 
in shocked regions. This exploration can be extended and made more realistic (in terms 
of upstream structure) by investigating internal evolution of shocked layers within larger 
fully-turbulent clouds having a range of mean Mach number; we are currently pursuing a 
numerical study along these lines. The models presented here may be thought of as inves- 
tigating self-gravitating structure growth within the first strong shocks to develop inside a 
cloud. 



^ Other recent simulations of post-s hock structure formation in converging; fl ows have similarly assumed 
uniform density for the inflow (see e.g. iHeitsch et al.ll2008t iBaneriee et al.ll2009L and references therein). 
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4. Development of structure and core evolution 



As discussed in Section 1, iGong fc Ostrikerl (120091 ) proposed a unified model for core 



formation and evolution in supersonic turbulent environments. Based on spherical-symmetry 
numerical simulations, four stages were identified: core building, core collapse, envelope infall 
and late accretion. The duration of each stage, and the structure and kinematics of cores 
at varying stages were also analyzed. While the comparison of those results to observations 
is very encouraging, the assumption of spherical symmetry is clearly unrealistic. One of 
the key goals of this work is to check if core building and collapse still develop in a similar 
manner when the spherical-symmetry assumption is relaxed. Because the time step becomes 
very short in late stages, we halt the simulations; thus the current models do not address 
envelope infall and late accretion stages. 

Figure [1] shows evolution of the surface density (eq. [36]) for models with = 1.1 (left 
column), Ai = 5 (middle column) and = 8 (right column), all with same realization for 
the perturbation velocities. The top panel of each column shows the surface density very 
early on; the patterns are identical but the amplitudes are different. The bottom panel 
shows the surface density when the most evolved core collapses for each case. Hereafter 
we shall use tcou to denote the total time to reach collapse of the most evolved core, in 
terms of the code unit to (eq. E]). The four images from top to bottom in the same column 
show the surface density at four instants: t = 0.001 to, l/3^coib 2/3i(:coii, and tcoU- Note 
that tcou = 0. 636^0, 0.280to and 0.232to for the Ai = 1.1,5 and 8, respectively. These three 
simulations have low initial perturbation amplitude (cf. eq. 



From Figure [U two features are immediately apparent. First, the input perturbation 
field patterns determine the later structural evolution and there is a "family resemblance" 
for the models at different Mach number. This is because the post-shock dense layer retains 
a memory of the perturbation velocity fields in the direction parallel to the plane of the layer 
since and Vy are unchanged across the shock interface. Comparing the first plot to the 
last plot of each column, cores form in regions where the density perturbation amplitudes 
are initially higher than the surroundings as a result of convergence in the x — y plane. 
These overdense regions develop into long, thin filaments, within which cores grow and then 
collapse. 

Second, the specific properties of cores, such as the total number and individual volumes 
(as well as their masses), are determined by M.. The dense cores for = 1.1 are smoother 
than the cores for = 8, and they cover larger areas. During the middle and late stages 
of evolution, more small scale filamentary structures are evident in the higher Mach number 
cases. At a given scale, the input and Vy perturbations are higher for larger Al, with 
the resulting compressions making more prominent "burrs" around cores. The "burrs" are 
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also less smoothed for the high Mach number cases, because the shorter free-fall time at the 
higher post-shock density means that the core collapses sooner. Thus, as the velocity of the 
converging flow and additional perturbations increases, the result is smaller, denser, more 
irregular, and more "hairy" cores. 

Figure [2] shows evolution of surface density and the mean in-plane velocities {vx) and 
{vy) for the M. = 'b model shown in Figure [H The mean velocities are calculated by (f ) = 
J pvdz/ J pdz with v = Vx or Vy. The left column shows surface density, and the middle 
and the right columns show {v^), (vy) respectively. At early stages, only scattered high 
surface density spots appear. The large-scale spatial correlation of these overdense regions 
is evident, however, even at early times. The mean velocities also have small amplitudes at 
early stages. The large-scale converging (in-plane) velocity regions that eventually lead to the 
most prominent filaments are already evident from the first frames, however. At late stages, 
the overdense regions start to collect into filaments. The converging (in-plane) velocities grow 
due to self-gravity of the forming filaments; in a ddition, purely hydrodynamic insta bilities 



(such as the nonlinear thin-shell instability, e.g. IVishniadll994 



Heitsch et al.l 120071 1 in the 



shock-bounded layer may enhance early growth of perturbations]^ When converging in-plane 
flows become supersonic, discontinuities in the density and velocity develop. These sharp 
fronts, as well as the collapsing motions centered on the most evolved cores, are evident in 
Fig. |2]at t = 11/12W,W. 

Thus, we see that turbulent motions even at sub-pc scales seed the growth of structures, 
and self-gravity reinforces and amplifies these motions. The growth of dense cores and larger 
scale filaments is simultaneous, both a consequence of turbulence and self-gravity. 



Similar to our results in iGong fc Ostrikerl ( l2009l ) for spherical symmetry, we find that 
core building lasts most of the time up to tcou, while the core collapse itself is rapid for the 
most evolved cores. Defining the "supercritical" period as the stage at which Pcenter/Pedge > 
10 for the most evolved core, this first occurs at 0.589^0,0.240^0 and 0.209 to respectively 
for the Ai = 1.1, 5 and 8 models shown in Figured] (we note that pedge is close to the post- 
shock density). Taking the difference with tcoib ^^supcrit/^o = 0.047,0.040 and 0.023. From 



Gong &: Ostrikerl (120091 ) . the supercritical stage lasts about 10% of tcoii for cores found in 
shocked converging spherical flows. For the three cases shown here, Atsupcrit /^coii is 7%, 14%, 
and 10%, consistent with our previous results. The core building stage lasts about 90% of 



To express At 



supcrit 



in terms of observables, we renormalize using the mean core den- 



We have conducted comparison tests of selected models without self-gravity, finding that surface density 
fluctuations can grow to order-unity level in high Mach number cases. 
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sity pmcan at the instant of collapse. This quantity, Atsupcrit/i//(Pmcan) = Atsupcrit/^o x 
3.27(pmean/Po)^''^ IS measured to be 0.9,2.1 and 0.8 for Ai = 1.1,5 and 8 respectively; i.e. 
Atsupcrit is comparable to t//(pmcan)- The values of Atgupcnt are 6.6 x 10^ yr,5.6 x 10^ yr 
and 3.2 x 10^ yr for Ai = 1.1, 5 and 8 respectively, if we take the inflowing ambient medium 
density as n^.o = 100 cm~'^; these are reduced to 2 x 10^ yr, 1.7 x 10^ yr and 1 x 10^ yr for 
HHfi = 1000 cm~^. 

Figure [3] shows the cross-sections of the density and velocity field across the center of 
the most evolved cores (the locations of these cores are indicated in Figured]) for = 1.1,5 
during the late collapse phase. The instants of the plot for J\4 = 1.1,5 are 0.625 to and 
0.273 to respectively. The top panels show the x — y cross-section of density and velocity 
vectors composed of and Vy in the same plane. The bottom part shows the x — z cross- 
section and velocity vectors composed of and Vz- The velocity field clearly shows inward 
collapse. The amplitudes of the velocity field are smaller in the outer part and larger in the 
inner part, indicating the core is at a very late stage of the "outside-in" collapse. 

Figure m and Figure E] show the evolution of the density and velocity profiles of the cores 
in Figure O The density profiles are azimuthally-averaged over the x — y plane. The velocity 
profiles are along each cardinal axis (x, y, z) through the core center. The instants for the 
four profiles have equal intervals 0.027to for = 1.1 and equal intervals 0.019 to for = 5 
respectively. The first instant for both cases is subcritical (i.e. Pcenter/Pedge < 10) and the 
second instant is close to tsupcrit- The dramatic increase of the central density during collapse 
is clearly evident for both cases, and the collapse develops in an "outside-in" manner with 
the maximum in v moving inward in time. The density profile approaches the asymptotic 
"Larson-Penston" profile p/po = 8.86(r/Lj)^^/(27r)^ at the instant of central singularity for- 
mation, and the in-plane velocities Vx-,Vy approach — 3.3cs, which is the "Larson-Penston" 
limit. Before the time tsupcrit is reached, the velocity is subsonic throughout the core re- 
gion. For all of the simulations we have conducted, the peak of the velocity profile becomes 
supersonic only at the very end of the collapse stage, similar to the results shown here. 

Overall, we conclude that the evolution of individual cores in these 3D simulations 



follow s a similar progression to the spherically-symmetric ID simulations of lGong &: Ostriker 



(120091 ). The core building stage lasts over 90% of the time to collapse, and cores become 
more stratified over time. The onset of the collapse is in an "outside-in" manner, and leads to 
a dramatic increase in the central density. As a central singularity is approached, the density 
and velocity profiles approach the "Larson-Penston" asymptotic solution. These cores form 
and collapse within larger-scale filaments that also grow in contrast over time. 
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5. Core-finding method 

The algorithm adopted for core- finding can e ither subtly or more seriously affect the 



core properties that result (e.g. iPineda et al.ll2009[ ). The most commonly-used methods in 



observational work are based on contouring colum n density or emission intensity (e.g. the 
popular Clump find method of IWilliams et al. I ll994h . For theoretical work, density-contouring 



methods, sometimes incorporating furthe r tests to determine i f a structure is gravitationally 



bound, have frequently been used (e.g. iGammie et al.l l2003l ). Here we shall instead use 



the gr avitational potential isosurfaces to identify cores. In very recent work, ISmith et al. 



( I2OO9I ) took a similar approach, noting that one advantage of the gravitational potential is 
that it yields smoother core boundaries than the density. Another advantage is that the 
gravitational potential connects more directly to the fundamental physics that determines 
core evolution. During formation stages, self-gravity gathers material to build up cores, and 
later it drives the collapse of supercritical cores. 

To identify cores via the gravitational potential, we first find and mark all the local 
minima of the gravitational potential; second, we find the largest closed potential contour 
(or isosurface) surrounding each individual minimum. In the second step, we increase the 
contour level from the bottom of a given potential well step by step until it violates another 
minimum's marked territory. We define the region enclosed by the largest closed contour as 
a core. The contour interval A$ has negligible effect on the results as long as it is small 
enough (typically < 0.03c^). If the distance between two minima is smaller than 10 pixels 
(corresponding to a physical distance ~ 0.03 — O.lpc for nn^ ~ 10^ — lO^cm"^), the regions 
associated with these two minima are merged and treated as a single core. Since we do not 
continue the simulation after the most evolved core collapses, we apply the algorithm to the 
last output from each simulation. 

Since gas with sufficient thermal and kinetic energy need not be permanently (or even 
temporarily) bound to a given core, the gravitational potential is not the final word. The 
lower density outer parts of a core are the most subject to loss. We can test this effect on 
core identification by adding thermal energy to the gravitational energy, and only assigning 
a given fiuid element to a core if Eth + Eg < 0. For any fiuid element, the specific thermal 
energy is taken to be Eth = 3/2Cg, and the specific gravitational potential energy is taken to 
be Eg = $ — $max5 where $max is the potential of the largest closed contour that defines the 
core. □ Including a thermal energy condition in core definition decreases the volume (or area 



''We note that \Eg \ for a core embedded within a dense filament (or sheet) may be much lower than \Eg \ for 
the same core in isolation. In assessing whether a core is bound, it is crucial to take tidal gravity effects into 
account. If these tidal effects are neglected, \Eg\ will be overestimated by a factor ^ I]corc/(Scorc ~ SfHamcnt), 
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in 2D) of the cores. Of course, the thermal energy can in fact be radiated away, so that gas 
that is initially near the largest closed contour may become more strongly bound after the 
interior of a core collapses. In this case, the potential alone could determine the final core 
mass. Short of following cores through the final stages of star formation, we consider it useful 
to compare cores with and without a thermal - gravitational energy criterion. Hereafter, we 
term our core-finding method "gravitational identification" (GRID). We refer to the region 
within the largest closed gravitational potential isosurface surrounding each local minimum 
as a GRID-core. For each GRID-core, the region which has E'th + -E^ < is referred to as a 
bound GRID-core. 

Because volume density data cubes are not directly accessible in observations, three- 
dimensional gravitational potential contouring is only applicable to model data from numer- 
ical simulations. It is therefore interesting to explore gravitational potential contouring of 
surface density maps, which are direct observables. To identify cores in a surface density 
map, we have to calculate the gravitational potential first. For a layer of half-thickness 
H, the gravitational potential component $fc, 20 of surface density component Sj, (Fourier 
transform of equation (I5B]) ) in phase space is 

*^'^'^ = -|fc| {i + \kH\y ^^'^ 

where \k\ = ^/kl'+k^■ Note that for \kH\ ^ 1, ^k,2D ~ —4:nGpk/k'^, which is the solution 
of the Poisson equation in three dimensions, for pk = Sfc/2if. For \kH\ ^ 1, eq. fH5l) is the 
solution of the Poisson equation for an infinitesimally thin layer. The gravitational potential 
$20(3;, y) is the inverse Fourier transform of $fc,2D- Given the 2D gravitational potential field 
^2D{x,y), we can apply the GRID procedure as for 3D. In Section 6, we will compare the 
results from GRID using ^{x,y,z) and $20(2;,?/) (using H = 6z). Hereafter we use "2D" to 
denote the results from applying the GRID method to surface density and "3D" for applying 
the GRID method to the volume density. 

As an example. Figure E] shows the comparison of GRID-cores and bound GRID-cores 
between 3D and 2D for Ai = 5 and 9. The top portion shows core areas identified for the 
J\4 = 5 model using <l> (top left) and $20 (top right). The bottom portion shows the same 
comparison for = 9 with cores found from $ (bottom left) and from $20 (bottom right). 
(Note that the Ai = 5 and Ai = 9 simulations have the same initial velocity perturbations 
patterns, which is why the overall structure is similar). In all plots, the areas enclosed by 
yellow contours are the GRID-cores and the areas enclosed by red contours are the bound 
GRID-cores. The core areas for the 3D plots are the projection of the 3D core volume onto 



which is quite large if the contrast between a core and its surroundings is modest. 
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the z = plane. For the M. = 5 model, the 2D and 3D core-finding procedures identify 12 
and 13 cores respectively; the cores and the bound regions are located at nearly the same 
positions. For the = 9 model, 7 cores are identified for both cases. One bound core in 
2D lacks a 3D counterpart, implying the corresponding potential well in 3D is too shallow 
(see discussion of potential well depths in Section 6). 

In addition to finding almost all of the same core centers (defined by the potential 
minimum), the areas marked by the 3D and 2D GRID algorithms are almost the same. 
Figure [7] show the results of GRID for four simulations for Ai = 5. The white contours 
mark GRID-cores from 3D density and the green contours mark GRID-cores from 2D surface 
density. The red and yellow contours mark the bound GRID-cores for 3D and 2D respectively. 
The areas identified for the cores agree quite well. Over all, we conclude that the 2D GRID 
algorithm can give nearly identical core-finding dbTQSbS clS the 3D GRID algorithm. 

In spite of the overall similarity between 2D and 3D GRID-core finding, there are minor 
differences in the results. In the each panel of Figure [TJ a few GRID-cores in relatively low 
density regions are identified in 2D but not in 3D. In comparing core properties between 2D 
and 3D, we shall apply additional resolution criteria to eliminate these small, shallow cores. 

6. Core properties 

To obtain a sufficient statistical sample, we conduct 20 simulations for each value of 
the Mach number {M. = 1.1, 2, 3, 4, 5, 6, 7, 8, 9) and compute GRID-core masses and radii for 
each model (180 models total). Each of the 20 simulations for a given Ai is perturbed by 
a different realization of the velocity field. As an example of the differences with different 
random realizations of the power spectrum. Figure [7] shows the snapshots of surface density 
at a late stage for four different Ai = 5 simulations. The 3D GRID core numbers are 9, 
6, 9 and 7. The corresponding core mass ranges are [0.00151,0.158] Mo, [0.0051,0.128] Mq, 
[0.0013, 0.242] Mq and [0.031, 0.250] Mq. The core numbers and core masses from simulations 
with different seeds are in a similar range; the same is true for cases with other Mach numbers. 

The GRID-core masses for 3D and 2D are M^d = j p dxdydz and M2D = / ^ dxdy, 
respectively. The GRID-core radius for 3D is defined as the equivalent radius of a 3D sphere 
with the same volume Vsd: ^"30 = (3V3D/47r)^/'^. The effective 2D GRID-core radius is 
calculated from the area 5*20 of the core region as: r 2d = {820/ '^Y^'^- To ensure that 
identified GRID-cores are numerically well-resolved, we only retain cores with effective radii 
> 4 zones. We define a background surface density as the mean of the bottom 10% of the 
surface density; this mean value can be subtracted from the surface density in the core region 
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when calculating M2d- As mentioned in Section 2, a more restrictive definition includes 
only gas with thermal plus gravitational energy negative; these bound GRID-cores are first 
identified by the gravitational potential, and then pixels are excluded if the sum of thermal 
energy and gravitational potential is greater than 0. 

Figure [8] shows M2D versus M3D for GRID-cores, for each Mach number of the low 
amplitude perturbation set. Note that only cores with same center of the local potential 
minima are shown here. Both 2D GRID-core masses without background subtraction (M2D, 
diamonds in the figure) and 2D GRID-core masses with background subtraction (M2D,bs5 
dots in the figure) are shown versus M3D. For large masses, M2D agrees well with M3D while 
-^2D,bs is slightly lower than M3D. For small masses, M2D,bs agrees better than M2D with 
M3D. Both M2D and M2D,bs agree with M3D better for high mass than low mass. 

Figure M shows a similar comparison of bound GRID-cores for 2D and 3D. The back- 
ground surface density is subtracted for 2D GRID-core masses, so that we show M2D,bs,th 
versus M3D,th- Here, the subscript "th" represents inclusion of a thermal energy criterion 
in defining bound GRID-cores, which eliminates most of the small cores. At high masses, 
^2D,bs,th agrees with M3D,th for bound GRID-cores better than M2D,bs agrees with M3D for 
the whole set of GRID-cores. This is because only zones sufficiently near the potential mini- 
mum where Eth+Eg < are included in bound GRID-cores; these regions are not sensitive to 
projection effects. At low masses, M2D,bs,th exceeds M3D,th for bound GRID-cores, meaning 
that imposing the thermal - gravitational energy criterion affects M3E) th niore than M2D,bs,th- 

To understand the difference between the 2D and 3D GRID-core masses, we consider 
the shape of the gravitational potential well for surface density and volume density. From 
equation (H5|) . $2D,fc oc —k~^ whereas $3d,a; oc —k~^. At larger fc, corresponding to smaller 
scales, |$3d| decreases faster than |$2d|- That means that the small 2D GRID-cores cover 
more area than small 3D GRID-cores, evident at the low end of each panel in Fig. [HI 
Also, gravitational potential wells of middle-sized 2D GRID-cores are deeper than those of 
3D middle-sized GRID-cores. If the shallow parts of the potential are excluded by apply- 
ing a thermal energy requirement, 3D GRID-cores are affected more than 2D GRID-cores. 
Moderate-mass GRID-cores that have M2D,bs and M3D comparable will thus have M3D,th 
lower than M2D,bs,th; as is evident in Fig. O As mentioned in Section 5, we include the term 
\k\H to allow for the non-zero thickness of the layer perpendicular to the plane. This can, in 
principle, help decrease the gap between the 2D and 3D gravitational potentials. In practice, 
however, we find that the value for H to make the central-to-edge value of $20 comparable 
to that for $ is smaller than 6z. Although the 2D and 3D gravitational potentials are not 
exactly the same. Figure [9] shows that 2D and 3D bound GRID-cores masses are generally 
close down to ~ IO^^Mq (which is < 1 M© for typical conditions, from eq. [5]). 
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Figure [To] shows histograms for the distributions of M2D,bs and M3D (all GRID-cores) for 
each A^, while Figure [TT] shows the histograms of M2D,bs,th and Mso.th (bound GRID-cores), 
both for low perturbation amplitudes. The distributions of M2D,bs and M3D are quite similar 
for all Ai, except slightly more low mass cores are identified for 2D at large Ai. When the 
thermal - gravitational energy condition is included in defining cores, the low-mass end of 
the distribution is removed; in Fig. [TTl the 2D bound GRID-cores have almost exactly the 
same distributions as 3D bound GRID-cores. 

Figure [T2] (all GRID-cores) and Figure [T3] (bound GRID-cores) show the median core 
mass (squares in figures) versus Ai from Figure [10] and [111 respectively. (We do not measure 
the peak because some of the histograms are irregular.) Figure [H] (all GRID-cores) and 
Figure [TS] (bound GRID-cores) show the same median mass - Ai relation for high amplitude 
initial perturbations. The breadth of the distributions at each Ai is indicated by vertical 
bars: the lower bar is the difference between the median and the first quartile, and the higher 
bar is the difference between the third quartile and the median. In Fig. [121 [IS] and Fig. 
[HI [151 we overlay hues showing the predicted critical mass at late stages (eq. [22] or [23] 
dashed line with M oc Ai~^), and the prediction for the mass that has grown the most at 
early time (eq. [2S]or[nil dot-dashed with M cx A^~^/^). The post-shock Bonnor-Ebert mass 
(M oc A^~^ from eq. [7]) is similar to the late-stage critical mass. 

As the Mach number increases, the post-shock density p ~ poAi"^ is higher. This lowers 
the Jeans length (as well as the Jeans mass and Bonnor-Ebert mass), permitting smaller 
(but denser) cores to form at high Ai compared to low Ai. However, high mass cores can 
still form at high Ai, as is evident in Figure [TO] and [TT] and the quartiles shown in Figures 
[T2] - [T5] at high At, the histograms extend to low mass, but the high mass part of the 
distribution is still present. This is consistent with the expectation that any scale above the 
critical scale can grow more nonlinear due to self-gravity (see eqs. [T^- [TT]) . 

Based on Figures [12] - [TS] we also note that the median mass versus Ai relations are 
quite similar whether cores are identified with the 2D or 3D gravitational potential. This 
is true for low or high amplitude perturbations, for both all GRID-cores and bound GRID- 
cores. This evidently shows that 2D cores have similar statistical properties to the 3D cores. 
Since the GRID algorithm is easy to implement for observational data, it appears to be a 
promising method for finding coresE 

Median masses for GRID-cores decline with increasing Mach number for both low and 
high amplitude perturbations (see Figs. [I21[I1|). These median masses generally lie above the 



'"^An IDL implementation of our GRID-core algorithm for use with observed data (FITS files containing 
surface density maps) is available from the authors. 
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values predicted from equations ([7]), f|T9|l and fl22l) (M oc A^"^) at late stages and below the 
values predicted from equation (|28|) (M oc A^^^/^) at early stages. The median GRID-core 
masses for high amplitude perturbations are slightly smaller than those for low amplitude 
perturbations, and the range of core masses for a given Mach number are larger. This reflects 
the fact that the percentage of small cores is higher when the perturbation amplitudes are 
higher. GRID-cores are identified based on the gravitational potential, and this potential 
reflects density structure, which arises from both turbulent and gravitational processes. Even 
without gravity, smaller scale masses would be expected in the higher- models because 
of their high turbulent amplitudes. For our simulations, the input perturbation amplitude 
at scale / is (5fiD(/) = (//Lj)^''^(A^/3)^/^ at 100% amplitude of perturbation (cf. eq. 
H3|) . Structures at scales I for which turbulent perturbations are supersonic will, even in the 
absence of gravity, be more prominent than those at smaller scale. For our adopted scaling 
of input perturbations with J\4, the sonic scale varies as /sonic oc Lj/J^, so that the mass at 
the sonic scale varies oc E(t) /g^j^;^. With S(t) oc A^tcou and tcou oc A^~^/^ (see eq. 129) and 
below), this predicts Mgonic oc A^~^/^. For later time t ~ tj (comparable to the flow crossing 
time for a cloud with a^n = 1^2), Msonic oc A4^^. Thus, the sonic mass scale, and hence 
the mass scale of nonlinear structures induced purely by turbulence, is expected to decline 
with increasing A4. 

For bound GRID-cores, the median mass vs. Ai decreases and then increases, for low 
amplitude perturbations (Fig. [13]), and is nearly flat for high amplitude perturbations (Fig. 
[T5]) . The high median mass at high Ai for bound GRID-cores may be due to a combination 
of effects, including numerical resolution and nonlinearity. The characteristic scale for self- 
gravitating perturbations decreases with increasing Mach number (either as r oc Ai^^^"^ 
for the most-grown core or r oc Ai~^ for critical perturbations; see Section 2). At high 
Ai, this may approach or fall below the minimum scale rmin = 4 zones = 0.016Lj that 
we require for the GRID-core radius to be well resolved. Since the post-shock density is 
oc A^^, the GRID-core mass would then increase at least oc At^r'^^^ at sufficiently high At. 
In addition, larger-scale, higher-mass regions initially have higher amplitude perturbations 
than smaller-scale regions, because of the input power spectrum with 5v oc Z^/^. If this 
initial "head start" allows the larger, more massive cores to become highly nonlinear before 
more rapidly-growing smaller-scale cores, the more massive cores will collapse (halting the 
simulation) before the lower-mass cores become strongly concentrated (with Ett < \Eg\) 
internally. With implementation of sink particles such that the simulations need not to be 
halted when the most evolved core collapses, and \Eg\ can grow for low-mass cores, it will 
be possible to test whether the median mass of bound cores decreases with increasing Ai, 
similar to Figs. [T^] and HH 

Figure [16] shows the GRID-core radii (as defined in Section 3) versus Mach number. 
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and Figure [T7| shows the bound GRID-core radii versus Mach number; these are for cases 
with low amphtude initial perturbations. Overall, the median radii for all GRID-cores and 
bound GRID-cores decrease towards higher Ai. This is consistent with expectations: high 
Mach number yields high post-shock density, and hence a smaller Jeans length; in addition, 
the higher amplitude of input turbulence at higher Ai makes the sonic scale smaller. The 
prediction for core radius based on turbulence alone would be the sonic scale from Equation 
( H3|l : Teff oc /sonic oc Lj/ Ai. The first core to collapse is predicted to have Xm oc A^~^/^ 
from equation ( l30ll . For late-time fragmentation, the relevant scale is the Jeans length in 
post-shock gas, which varies oc A^~^. For GRID-cores, the slopes are between these values, 
equal to —0.95 ± 0.13 for rcff,2D,bs and —0.72 ± 0.07 for rgfr^sD, for low amplitude initial 
perturbations. For bound GRID-cores, the power-law fit for median radius as a function of 
Mach number gives slope — 0.67±0.10 and —0.61 ±0.08 for 2D and 3D respectively. These are 
comparable to the result oc Ai~^^'^ from Equation f l5U]) . Although the overall slopes are 
close to —0.5, we note that the relation flattens at > 5, possibly due to our requirement 
that the effective radius must exceed 4 zones, or because the initial power spectrum favors 
larger cores. 

Figure [TS] shows the median collapse time of the most evolved core vs. Mach num- 
ber, for both low and high amplitude initial perturbations. They both follow power laws 
close to tcou oc Ai~^^'^, consistent with the time scale (see eq. [29]) predicted for growth of 
self-gravitating modes up to a given amplification Fmax- The coefficients for low amplitude 
initial perturbations and high amplitude initial perturbations are 0.69 and 0.51, respectively, 
compared to 0.34 from equation fl2^ taking Fmax = 1- With high amplitude initial pertur- 
bations, cores collapse earlier because the seed perturbations need not grow as much. Note 
that the naive expectation based on the Jeans time, taking Ppost-shock oc A^~^, would yield 
a steeper dependence t oc Poost-shock oc A^~^. Based on Fig. [181 it is evident that the first 



-shock 

cores in higher Ai cases collapse when the layer as a whole is only barely self-gravitating 
(^couAo ~ 0.2 — 0.3, compared to t^g ~ 0.22to from eq. [TT]), whereas the layer is more strongly 
self-gravitating at the first collapse for low-A^ cases. 

The shape of a core can be characterized by the eigenvalues of the moment of inertia 



tensor Jjj = J pXiXjd?^ (e.g. iGammie et al.ll2003l : iNakamura fc Lill2008[ ). Let a, 6 and c be 



the lengths of the principal axes and a > h > c. Then a prolate core has h/a = c/a, and 
an oblate core has h/a = 1. We have computed the moment of inertia and aspect ratios 
for all the cores identified in our simulations. For example, the aspect ratios of the most 
evolved cores shown in Figures [T] and [2] are h/a = 0.39, c/a = 0.25 for the Ai = 1.1 model 
and h/a = 0.28, c/a = 0.25 for t he Ai = 5 model. T hey are both (approximately) prolate 



according to the classification of iGammie et al.l ( 120031 ) 
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Figure [19] and Figure [20] show the distribution of core aspect ratios for each Ai for low 
and high amphtude initial perturbations respectively. Open circles represent GRID-cores, 
and dots represent bound GRID-cores. These distributions show a number of interesting 
features and trends. First, only a small portion of cores are oblate for each Ai, for both low 
and high amplitude perturbations. Second, more oblate-like cores appear when the initial 
perturbation amplitudes are higher. For low amplitude perturbations, at Ai = 1.1 and 2, 
c/a and b/a are mostly < 0.5, i.e. approximately prolate. But at larger M. for low amplitude 
initial perturbations, and all A4 for high amplitude perturbations, there are many cores in the 
triaxial and oblate regions. Also, large and massive cores tend to be more prolate. For low 
amplitude perturbations, at Ai = 1.1, almost all the cores formed are prolate and no small 
cores form (compared to high Mach number cases). The reason that the distribution is more 
oblate for higher amplitude perturbation (large Ai for low amplitude initial perturbations, 
and all Ai for high amplitude initial perturbations) is that more of the cores are at earlier 
stages of evolution. Figure [1] shows development of cores for Ai = 1.1,5 and 8. As is 
particularly clear for the stages shown in the = 1.1 model, structures are more oblate 
during the core-building stage than during the collapse stage. Cores evolve to become prolate 
when they collapse because the collapse happens first in the directions perpendicular to the 
larger scale filaments. For Ai = 1.1,2 models with low amplitude perturbations, only large 
cores form and they have evolved to the collapse stage and become prolate. Models with 
higher amplitude perturbations have a greater percentage of small cores that have not yet 
collapsed. 

We can also examine the relationship between core structure and kinematics in our 
simulations. Figure [21] shows the projected density field, velocity field and the velocity 
dispersion field along the line-of-sight for the Ai = 5 model shown in Fig. [6] We "view" the 
simulation at angles 0°, 30° and 60° with respect to the z axis, tilting toward the x-axis. The 
white contours mark the regions identified as GRID-cores, and the orange contours mark 
the bound GRID-cores. The projected density field is smeared as the tilt angle 9 increases. 
Since (fios) = (vx) sin(6') + {vz) cos{9), with {vz) = and the contribution from {vx) small at 
6 small, no obvious pattern is seen for (fios) at = 0° and 30°. At = 60°, when the (vx) 
contribution becomes larger, converging flow patterns similar to those seen in Fig. [2] become 
apparent, especially surrounding the diagonal line of small cores. As previously discussed, 
converging flows in the x-y plane create this high density filament, which then fragments 
into small cores. 

As Figure [2T] shows, the dispersions of the line-of-sight velocity of high density regions 
are generally subsonic, and are even smaller in the cores. Velocity dispersions are low in high- 
density regions for two reasons. First, if filaments lie between supersonic converging flows 
in the x-y plane, then post-shock velocities within the filaments will be subsonic. Second, 
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weighting by density picks out regions that are physically small along the line-of-sight. The 
increase of linewidth with size means that if a region is smaller than its surroundings along 
the line-of-sight, then the linewidth will be smaller than that of its surroundings. Thus, from 
a combination of low post-shock velocities (in the x-y plane), and spatially-limited scale (in 
the z direction), aios is low in filaments and lower in cores, as seen in Fig. [211 



7. Summary and discussion 

Stars form in GMCs pervaded by supersonic turbulence, and core formation theory must 
take these supersonic turbulent flows into account. In this work, we explore the physics 
of core formation in a dynamic environment, focusing on post-shock layers generated by 
collisions of supersonic flows. The framework we adopt - three-dimensional planar converging 
flows containing multi-scale turbulence - enables us to analyze the internal structure and 
kinematics of cores, and to investigate the relation between core properties and the inflow 
Mach number M.. We consider a range Ai = 1.1-9, and conduct 180 simulations with 
different realizations of the initial turbulent power spectrum, in order to obtain a sizable 
statistical sample. In addition to core masses and sizes, we measure aspect ratios. To 
define cores, we introduce a new method based on the gravitational potential, and compare 
properties of cores identified using $ (from the volume density) and $20 (from the plane-of 
sky projected surface density). 

Unlike previous studies of core evolution that begin with pre-existing cores, the present 
models include formation stages. Our initial density is uniform everywhere, and cores grow, 
via self-gravity, from turbulence-induced perturbations within the post-shock layer; when 
the Mach number is high, initial growth of density perturbations is aided by shock-driven 
hydro dynamic in s tabilit ies. Based on a set of spherically-symmetric numerical simulations. 
Gong Sz Ostriker proposed four stages for core evolution in dynamic environments: 



core building, core collapse, envelope i nfall, and late ac c retion . The key features during 



core building and collapse described in iGong &: Ostrikerl ( 120091 ) are verified here, for more 



realistic geometry. As the supersonic flows converge in a plane, two reversed shocks propagate 
outwards. With its high mean density, the stagnation layer between these two shock fronts 
becomes an incubator for self-gravitating cores. When these cores become sufficient stratified, 
they collapse. We halt the simulations at the instant of singularity formation in the most 
evolved core, because the time step becomes very short. 

Based on the analysis of our simulations, our chief conclusions are as follows: 



1. Cores with realistic properties are able to form in post-shock dense layers within 
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turbulent GMCs. Core building to become supercritical takes ~ 10 times as long as the sub- 
sequent "outside-in" collapse stage, which lasts a few x 10^ yr. The d uration of the supercrit- 
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2008; 


Evans et al. 


2009) 



2. At the time of singularity formation, the radial density profile within cores approaches 
the Larson-Penston asymptotic solution p = 8.86Cg/(47rGr^) and the velocity approaches 
the Larson-Penston limit — 3.28cs. This is consistent with previo us studies of spherical core 



collapse (see Section 1 for references). iTilley fc Pudritzl (120041) also found tha t p oc 



in their most massive cores, for turbulent simulations. As in iGong &: Ostrikerl ((20091), we 
therefore conclude that the Larson-Penston asymptotic solution is an "attractor" for core 
collapse, no matter how the collapse is initiated. 

3. Prior to collapse, the velocities within dense cores remain subsonic, in spite of the 
highly-supersonic flows that create them. This is true both for the ordered inflow, and for 
the mean internal velocity dispersion. This result is consistent with observations that most 



cores have subsonic non-thermal velocity dispersions ( MversI 



Caselh et al.ll2002l : iTafalla et al.ll2004j ; Kirk et al.ll2007t lAndre et al. 



1983 



Goodman et al. 



2007 



Lada et al. 



1998 



20081). 



The velocity dispersion can increase quite sharply at the edge of the core in ou r models (see 
Fig. [2T]) . intriguingly similar to a sharp transition seen in NH3 observations by lPineda et al. 
feoioh for the B5 core in Perseus. From some orientations, velocity dispersions in filaments 
containing cores may also be lower than in the surrounding gas (cf. Fig. I^T]) . 

4. At sub-pc scales, turbulent velocity perturbations (whether super- or subsonic) induce 
density perturbations that can grow strongly if the density is high enough for self-gravity 
to be important. In post-shock layers, turbulence and self-gravity collect gas into long, thin 
filamentary structures at the same time as the highest density regions within the filaments 
grow to become centrally-condensed cores. These filamentary structures containing embed- 
ded cores ar e similar to the structures in the Aquila rift and Polaris Flare clouds observed 
by Herschel jAndre et al.lboid iMen'shchikov et al.l[201ol ). 



5. Using the gravitational potential to identify cores is advantageous because it enables 
a core definition based on dynamical principles. For numerical simulations, the gravitational 
potential may be computed from the volume density (yielding $) or from the projected sur- 
face density (yielding $2d)- We show for our models that cores defined using $ and $2D are 
nearly the same, both for GRID-cores (defined by the largest closed potential isosurfaces) 
and bound GRID-cores (which additionally require i^^th + Eg < 0). Since $20 can be com- 
puted for observed clouds, using potential contours offers a promising new core identification 
method for application to high-resolution molecular cloud maps. IDL code implementing our 
GRID-core algorithm, suitable for application to observed data, is available from the authors. 
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6. We find that the range of core masses that forms increases as the Mach number 
Ai increases. Physically, this is because a larger range of spatial scales has significant 
perturbations when the turbulence amplitude is higher, and because the minim um mass to 



be gra vitationally unstable decreases as the density in the shocked layer increases. iBasu et al. 



( I2OO9I ) also found broader mass distributions when the turbulent amplitude is increased. At 
high Mach number, GRID-core masses range between ~ 10~^ - IMj, corresponding to ~ 0.05 
- 50 Mq for typical GMC conditions. 

7. Analytical arguments (see Section 2) suggest that the first core to collapse will have 
mass M oc A^~^/^, and that at late times, the minimum mass core will vary as M cx: Ai~^. 
Our numerical results for median core masses as a function of Ai lie between these two 
relations. When the core definition includes the condition that i?th + Eg < 0, the median 
mass increases at the largest Mach number. This may be due to the nonlinear "head start" 
of massive cores, such that lower mass cores have not yet become concentrated when the 
first core collapses (and the simulation is stopped). 

8. Analytical arguments (see Section 2) suggest that the effective core radius will decline 
with increasing Mach number, with powers between r^s oc Ai^^^"^ and r^s oc A^^^. Our 
numerical results show a decrease of with Ai in this range. For bound GRID-cores 
{Eth + Eg < 0), the relation is shallower than for GRID-cores defined by gravitational 
potential alone. 

9. The time for the first core to collapse in our simulations depends on Mach number, 
with tcou oc A^~^/^, and a slightly smaller coefficient for high-amplitude initial perturbations 
(see Fig. [18]). This scaling is consistent with analytic predictions for gravitational instability 
in a shocked converging flow (see eq. |29]). For high J\4, as is observed in GMCs, the first 
cores could collapse within a few Myr of cloud formation. For high Ai, the first cores collapse 
when the shocked layer containing them is only barely self-gravitating; this suggests that 
collections of stars can begin to form individually before they collapse together to create a 
cluster. 

10. A very small portion of cores are oblate, while most cores are prolate or triaxial. 
Large cores are preferentially prolate. The triaxiality of most c ores is consistent with previous 



resul t s from turbulent hydrodynamic and MHD simulations (IGammie et al.l 12003 



Li et al. 



2OO4J : iNakamura fc Lil l2008t lOffner et al.l 120081 ) . We also find that core shapes change as 



they evolve, from more oblate during early stages to more prolate during collapse. For high 
initial perturbation amplitudes, the distributions have a higher proportion of oblate cores 
because small cores are less evolved (at the time the first core collapses), compared to those 
in models with low initial perturbation amplitudes. 
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As noted above, the current models have provided evidence that the masses of cores 
that form depend not just on the mean Jeans mass in a cloud, but also on the cloud's 
level of internal turbulence at large scales, Equations fl22]) and fl25]) suggest that at 
late times, the characteristic core mass will follow oc ct^^Pq ^^^T^, where po is the mean 
density in the cloud. For the current simulations, however, we halt at the instant when the 
most evolved core collapses (because the time step becomes very short). This limits the 
condensation of small cores; they are present, but not yet strongly bound. In order to fully 
test the dependence of on cloud parameters, i t is necessary to implement sink particles 



[e.g. iKrumholz et al.l 12004 : iFederrath et al.l 120101 ) so that the simulation can run until all 
the "eligible" cores in the post-shock region have had the opportunity to collapse. Including 
sink particles, as well as studying shocked converging flows within larger turbulent clouds 
via mesh-refined simulations, represent important avenues for future research. 



We are grateful to Lee Mundy and Alyssa Goodman for stimulating conversations, and 
to the referee for a helpful report. This work was supported by grants NNX09AG04G and 
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Fig. 1. — Evolution of surface density projected in the z direction (color scale logS/So = 
logN/NQ, see eqs. [Ml ETj) for converging-flow Mach number = 1.1 (left column), = 5 
(middle column) and = 8 (right column) models with the same initial perturbation 
patterns. The four panels from top to bottom in the each column show surface density 
snapshots at four instants: t = 0.001 to, l/3tcoih 2/3tcoib and tcoih with tcou the duration of 
the whole simulation. These three simulations have 10% initial perturbation amplitude (see 
eq. HHl). The values of tcoii are 0.636to, 0.280to and 0.232to for A4 = 1.1, 5 and 8 respectively 
(see eq. [6] for definition of to)- Cores are clearly smaller and more irregular for high-Al 
models. The squares indicate the most evolved cores for = 1.1 and 5. 
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Fig. 2. — Evolution of surface density (left column, log color scale) and the in-plane velocity 
components (vx) (middle column) and (vy) (right column) projected in the z direction for the 
Mach number Ai = 5 model shown in Figured], where (v) = J pvdz/ J pdz. The four panels 
from top to bottom in the each column show four instants: t = 0.001 to? l/^Wi, ll/12tcoih 
and tcoib with tcou = 0.28to the duration of the simulation (see eq. |6] for definition of to)- 
In-plane velocity fields are initially low, but grow to become supersonic, creating filaments 
that fragment into cores. 
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Fig. 3. — Density and velocity field cross-sections at the time tcou in the most evolved core, 
for M. = 1.1 (left column) and Ai = 5 (right column). These correspond to the most 
evolved cores (as indicated with boxes) in Figured] for Ai = 1.1,5 respectively. The color 
scale represents x — y and x — z slices through the volume density (logp/po)- The direction 
and length of arrows indicate the direction and magnitude of the local velocity, with scale as 
indicated in the upper left. At this stage of collapse, velocities increase toward the center. 
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Fig. 4. — Radial density and velocity profiles during collapse, for the most evolved core 
shown in Figured] and Figure [3] for Ai = 1.1. The density profiles are averaged azimuthally 
in the x — y plane about the center of the core. The dashed line is the Larson-Penston 
asymptotic density profile p/po = 8.86(r/Lj)~^/(27r)^ (i.e. p = 8.86cl/[47rGr^]). The other 
three plots show the corresponding velocity profiles versus distance in the x, y and z direction, 
respectively. The instants shown are 0.549 to, 0.576 to, 0.603 to, 0.632 to ~ tcoii, with the most 
evolved profiles in each case having the largest excursions. The collapse develops in an 
"outside-in" manner with the maximum in v moving inward with time. The density profile 
approaches the Larson-Penston profile with time. 
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Fig. 5. — Same as in Figure H] for the most-evolved core of the A4 = 5 model shown in Fig. 
[U The profiles are shown at t = 0.219 to? 0.238 to? 0.257 to? 0.276 to? with the density at the 
final time reaching the Larson-Penston solution. 
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Fig. 6. — Comparison of GRID-cores defined via the gravitational potential computed from 
3D volume density ($, left column) and 2D projected surface density ($20, right column). 
The top row shows = 5 and bottom row = 9. The areas enclosed by yellow curves are 
the GRID-cores determined by the largest closed gravitational potential ($ or $20) contour 
surrounding a local potential minimum, and the areas enclosed by red curves are the bound 
GRID-cores. Color scale shows projected surface density (logE/Eo ) in all panels. Cores 
identified using $ and $20 agree quite well. 
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Fig. 7. — Late stage surface density (logS/So) and GRID-core comparison for four dif- 
ferent random perturbation realizations of the = 5 model. The snapshots are at 
t = 0.282to, 0. 304^0, 0.304to, 0.302to from left to right and top to bottom. The corresponding 
maximum densities are 1.0 x lO^po, 1.53 x 10^po;8-18 x lO^po, 1-34 x lO^po- The white and 
green curves are GRlD-cores defined by the largest closed contour of the gravitational po- 
tential ($ and $20 respectively) surrounding each potential minimum. The red and yellow 
curves are the bound GRID-cores obtained using $ and $2D; respectively. Except for a few 
small, shallow cores, the core-finding algorithms in 2D and 3D give quite similar results. 
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Fig. 8. — GRID-core mass obtained from from 2D (M2d) versus 3D (Man). Diamonds are 
M2D for 2D GRID-cores without background subtraction, and dots are M2D,bs for 2D GRID- 
cores with background subtraction. The mass unit Mq is given in equation ([5]). Sohd hues 
represent M2D = ^3^; higher-mass cores are consistent with this. 
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Fig. 9. — Bound GRID-core mass for 2D with background subtraction (M2D,bs,th), versus 
bound GRID-core mass for 3D (Af3D tii)- When the condition Eti^ + Eg < is included in the 
core definition, the lowest mass cores are eliminated and M2D,bs,th agrees well with MsD.th 
down to ~ 10~^Mo. 
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Fig. 10. — Histograms of all GRID-core masses found in all simulations for each Mach 
number M. for low amplitude perturbations. Solid lines are for 3D GRID-cores [Mzd] and 
dashed lines are for 2D GRID-cores with background subtraction (M2D,bs)- The 2D and 3D 
distributions are similar for all Mach numbers. 
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Fig. 11. — Same as in Figure [TOl except for bound GRID-cores (i.e mass is Ms^.th and 
^2D,bs,th)- When the condition E^y^ + Eg < is apphed, most of the low mass cores are 
ehminated, for every Mach number. The 2D bound GRID-cores have almost the same mass 
distribution as 3D bound GRID-cores. 
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Fig. 12. — Median GRID-core mass M versus Mach number M. of the inflow. The left 
panel is for 2D GRID-cores (M2D,bs) and the right panel for 3D GRID-cores (Msd). Vertical 
bars indicate quartiles of the distribution. Also shown is the expected mass dependence for 
early gravitational fragmentation given by equation (128|) (with M oc A^^^^^, dot-dashed), 
and late gravitational fragmentation given by equation ([22D (with M oc dashed). The 

critical Bonnor-Ebert mass at the post-shock density (see eq. [7j) is similar to the late-stage 
prediction (M oc A^^^, dashed). The relation between median core mass and M. is quite 
similar for 2D and 3D cores. Core mass declines with increasing Mach number A^, lying 
between the M oc (early stage) and M oc A^^"*^ (late stage) fragmentation predictions. 
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Fig. 13. — Same as in Figure [T^ but for bound GRID-cores [Et\i + Eg < 0, i.e. M is M2D,bs,th 
or M3D,th)- 
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Fig. 14. — Median GRID-core mass M2D,bs and M3D, as shown in Figure [121 but for 
high amphtude initial perturbations. The median masses are shghtly smaller than for low 
amplitude initial perturbations, but follow a similar trend. 



1.00 



0.10 



0.01 



2D 



3D 



10 1 



10 



Fig. 15.— Median bound GRID -core mass A^2D,bs,th and -M3D,th (i-e. i?th + Eg < 0) as in 
Figure [131 but for high amplitude initial perturbations. 
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Fig. 16. — Median GRID-core radius versus Mach number M. for low amplitude initial 
perturbations. Core sizes are defined using the largest closed contours of the gravita- 
tional potential in 2D ($2D; left) and 3D ($, right). Vertical bars indicate quartiles of 
the distribution. The dotted lines are power-law fits: refr,2D,bs/-^j = 0.23qj|A^~°'^^''=°'^^ and 
reff,3D/i^j = 0.16[]:}|A^- 
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Fig. 17. — Same as in Figure [16] but for bound GRID-cores (-Eth + Eg < Q). The power-law 
fits are reff,2D,bs,th/^J = 0.15[]i|A^-o-67±o-io and reff,3D,th/i^j 
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Fig. 18. — Collapse time of the most evolved core, tcoii, versus inflow Mach number M. for 
low amplitude (squares) and high amplitude (diamonds) initial perturbations. Each value is 
the median of tcoii for 20 simulations for each M.. Vertical bars indicate quartiles of these 
20 values of tcou- The solid line least-squares fits are: tcou/^o = 0.69A^^°''** (low amphtude) 
and tcoiiAo = 0.51A^^°'^^ (high amplitude). The scaling is comparable to tcou oc A^^°'^, as 
predicted by equation (1291) . The simulation time unit to, based on the mean inflow density, 
is given in equation ([6]). 
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Fig. 19. — Distribution of three-dimensional core aspect ratio for each Mach number for 
low amplitude initial perturbations. Cores lying on c/a = b/a are formally prolate and 
along b/a = 1 are formally oblate. We subdivide (see diagonal lines) and classify as follows: 
approximately prolate (between c/a = 1 and c/a = 1.56/a — 0.5), triaxial (between c/a = 
1.5b/a — 0.5 and c/a = 3b/a — 2) and approximately oblate (between c/a = 3b/a — 2 and 
b/a = 1). Open circles are GRID-cores defined by the gravitational potential contours alone. 
Dots are bound GRID-cores, with the additional requirement Eth + Eg < 0. 
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Fig. 21. — Observations of one of the = 5 models shown in Fig. |6] from different 
angles. The first column shows the surface density (color scale logS/Eo); the second column 
shows the line-of-sight velocity and the third column shows the dispersion of the line-of-sight 
velocity (linear color scale, in units of c^). The three rows from top to bottom show the 
observed fields for ^^tiit = 0°, 30° and 60° respectively. The white curves are the GRID-cores, 
and the orange curves are the bound GRID-cores. Note that core regions have low internal 
velocity dispersions. 



